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ABSTRACT 

A Bayesian statistical formalism is developed to quantify the level at which the mass-function 
(dN/dm cc m^") and the projected cumulative mass fraction (/) of (CDM) substructure in 
strong gravitational-lens galaxies, with arcs or Einstein rings, can be recovered as function of 
the lens-survey parameters and the detection threshold of the substructure mass. The method is 
applied to different sets of mock data to explore a range of observational limits: (i) the number 
of lens galaxies in the survey, (ii) the mass threshold, Miow, for the detection of substructures 
and (iii) the uncertainty of the measured substructure masses. We explore two different priors 
on the mass function slope: a uniform prior and a Gaussian prior with a - 1.90 ±0.1. With 
a substructure detection threshold M|ow = 3 x lO^M©, the number of lenses available now 
(«/ - 30), a true dark-matter mass fraction in (CDM) substructure < 1.0% and a prior of 
a = 1.90 ± 0.1, we find that the upper limit of / can be constrained down to a level < 1.0% 
(95% CL). In the case of a uniform prior the complete substructure mass distribution (i.e. mass 
fraction and slope) can only be characterized in a number of favourable cases with a large 
number of detected substructures. This can be achieved by an increase of the resolution and 
the signal-to-noise ratio of the lensed images. In the case of a Gaussian prior on a, instead, it is 
always possible to set stringent constraints on both parameters. We also find that lowering the 
detection threshold has the largest impact on the ability to recover a, because of the (expected) 
steep mass-function slope. In the future, thanks to new surveys with telescopes, such as SKA, 
LSST and JDEM and follow-up telescopes with high-fidelity data, a significant increase in 
the number of known lenses (i.e. »10'*) will allow us to recover the satellite population in 
its completeness. For example, a sample of 200 lenses, equivalent in data-quality to the Sloan 
Lens ACS Survey and a detection threshold of IO^Mq, allows one to determine / - 0.5 ± 0. 1 % 
(68% CL) and a = 1.90 + 0.2 (68% CL). 

Key words: gravitational lensing — dark matter — galaxies: structure — galaxies: haloes — 
methods: statistics 



1 INTRODUCTION 

In the context of the cold-dark-matter paradigm, a significant 
number of substructures, with a steep mass function, is expected 
to populate the dark halo of galaxies. In galaxies as massive 
as the Milky Way, for example, of the o rder of 10"* substruc- 
tures are predicted inside the virial radius jPiemand et al.l [20081 ; 
ISprin gel et al."2008'). although only about 20 have been so far ob- 
served (Zucker etal. 2004; Willman et al. 2005; Belokurov et al. 



20061 ; lGrillmaij|20Q 6; Martin et al. 20061; ISakamoto & Hasegawa 
20od ; IZucker e t a^ 2006 a.h; .Beloku rov et al.l l2007l; llbata etal 
20071 ; llrwin eT al. 2007; Maiewski et al. 2007; IWalsh et al.ll2007l : 



Zucker et al l2007l ; lBelokurov et alj2008i) . A clear comparison be 



tween the simulated and the physical reality, however, is strongly 
hampered by the difficulty of directly observing substructures in 
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distant galaxies, as well as in the Local Group. 
While major improvements in the observations and numerical sim- 
ulations have not yet significantly alleviated the satellite crisis, new 
techniques have been proposed for the indirect and direct detec- 
tion of subhaloes that may have a high mass-to-light ratio. In our 
own Galaxy, CDM substructures can be, in p rinciple, id entified 
via th eir effect on stellar streams l llbata et all[2002 ; M ayer et al.l 
l2002h or via th e dark matter annihila t ion signal from their high- 
densi t y centres i Bergstrom et alJI 19991; Calcaneo-Roldan & Moord 



l200d ; Istoehr et al.ll2003l ; IColafrancesco et al.ll2006h ; gravitational 
lensing, on the other hand, allows for direct detection (measure- 
ment of the substructure gravitational signature) in the central re- 
gions of galaxies through flux-ratio ano malies and distortions of 
extended Einstein rings a nd arcs, (e.g. Mao & Schneideil l998l; 
Metcalf & MadatJ 1200 ll ; IPalal & Kochanekl |2002| ; iKoopman^ 



2OO3) . Interestingly enough, results based on flux-ratio anomalies 
reverse the satellite crisis with a recovered mass fraction in sub- 
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structure which seems to be larger than predicted by nu merical 
simulations (e.g. lMao et al.l2004lMacci6 & Mirand j2006l) . While 
this discrepancy might not be easily accommodated by an increase 
in the resolution of simulations, the correct interpretation of the 
fl ux-ratio anomalies is s t ill sub ject of discussion. 
In lVegetti & Koopmans! l f2009h , we introduced a new adaptive-grid 
method, based on a Bayesian analysis of the surface brightness dis- 
tribution of highly magnified Einstein rings and arcs, that allows 
the identification and precise quantification of substructure in sin- 
gle lens galaxies. This technique does not depend on the nature of 
dark matter, on the shape of the main galaxy halo, strongly on the 
density profile of the substructure, nor on the dynamical state of the 
system. It can be applied to local galaxies as well as to high redshift 
ones, as long as the lensed images are highly magnified, extended 
and have a high signal-to-noise ratio. Unlike flux ratio anomalies, 
extended optical images are little affected by differential scatter- 
ing in the radio or microlensing in the optical and X-ray. If sub- 
structures are located close to the lensed images, the method allows 
the determination of both their mass and position, although as the 
distance between the substructure and the Einstein ring increases, 
the mass model becomes more degen erate. F inally, thanks to its 
Bayesian framework, the method of Vegetti & Koopmans ( 2 0091) 
requires hardly any subjective intervention into the modelling and 
any ass umption can be objectively tested through the Bayesian Ev- 
idence ( lMacKavlll992l) . 

In this paper, we show how results from data sets with more than 
one lens system (i.e. the number of detections and their masses) can 
be combined to statistically constrain the fraction of dark matter in 
substructure and their mass function, as function of the survey pa- 
rameters and limits. The combination of multiple data sets becomes 
important when trying to constrain the slope of the substructure 
mass function. In a single lens potential, in fact, the number of de- 
tectable substructures can, in certain cases, be as small as zero or 
one. More than one lens is, therefore, required in order to improve 
statistics and properly sample the mass function. On the contrary 
an upper limit to the mass fraction in substructure can always be 
set. 

The paper is organized as follows. In Section 2, we outline the sta- 
tistical formalism of satellite detection; in Section 3, application of 
the method is discussed and in Section 4 conclusions are drawn. 



2 BAYESIAN INTERPRETATION OF SUBSTRUCTURE 
DETECTIONS 

In this section, we outline a Bayesian formalism that allows us to 
give a statistical interpretation to the detection of dark substructure 
in gravitational lens galaxies and to recover the properties of the 
substructure population. Thanks to Bayes' theorem, the likelihood 
of measuring a mass substructure, can be translated into probabil- 
ity density distributions for the substructure mass fraction / and 
their mass function dN/din oc m"" (when dN/din is normalised to 
unity, we refer to it as dP/dm), as function of the mass measure- 
ment errors and the model parameters. The latter include the min- 
imum and maximum mass of dark matter substructure, Mmin and 
Mmax, between which the mass fraction is defined, and the lowest 
and highest mass we can detect, M]„„ and M^igh. More precisely, 
Miow should not be interpreted as a hard detection threshold, but as 
a statistical limit above which we believe a detection to be signifi- 
cant, at the level set by the mass measurement error cr„, (i.e. with a 
signal-to-noise ratio of M|ow/o",„). 



2.1 Likelihood of the substructure measurements 

We derive an expression for the Likelihood of observing sub- 
structures for a given sample of «/ lenses. We assume the cumula- 
tive dark-matter mass fraction of substructure /(< R/R^i^), within a 
cylinder of projected radius R of the lens, to be the same for all lens 
galaxies and the detections of multiple substructures in one galaxy 
to be independent from one another. We also assume that the num- 
ber of substructures populating a given galactic halo fluctuates with 
a distribution which is Poissonian. We know that not all these sub- 
structures are observable, but only those with the right combination 
of mass and position. 

The Likelihood of measuring substructures, each of mass m,, 
in a single galaxy, in general can be expressed as the proba- 
bility density of having substructures in the considered lens, 
times the normalised probability density P (nij, R\ p,a) of actu- 
ally observing the mass m, within the projected radius R, where 
p = {Mn,i„, Afnmx, Mow, A^hieh} is a vector containing all the fixed 
model parameters introduced above, so that 



£(n„ m I a, f,p) 



IA{aJ,<Rr' 



Y\P{m„R\p,a) 



(1) 



where m contains all the substructure masses ot, and yu(a, /, < R) 
is the expectation value of the number of substructures in a generic 
aperture with dark matter mass Mdm(< R), which will be discussed 
in more details in the following section. The second factor in equa- 
tion (O describes the likelihood with which a substructure can be 
identified as a function of its mass m, and position R. Although the 
specific shape of P (ot,, R\ p,a) might change from one system to 
another, its overall trend is essentially the same for all lenses: high 
mass perturbations, on or close to the lensed images are the most 
likely to be confidently observed. For each considered lens system, 
P(inj,R I p,a) can be reconst ructed by Monte Carlo exp lorations 
of the model parameter space ( I Vegetti & Koopmanj|2009l) . 
Here, for the sake of simplicity, we assume that the probability of 
measuring a mass nij for a single substructure is a Gaussian func- 
tion independent of position. This is indeed true for small regions 
around the Einstein ring of the lensed images, so that P (m, , R\ p,a) 
reduces to 



P {nii \p,a) 



JM„;„ dm \i 



V2jrtr,„ 



-dm 



Jmi™., JiWmi,, dm I, 



-(m-m' )- llm- 



(2) 



dm dm' 



where the Gaussian convolution expresses the scatter in the sub- 
structure mass due to a measurement uncertainty and dP/dm\,r„e is 
the true substructure mass function as defined in equation (6). The 
assumption of Gauss ian errors on the substructu re mass is moti- 
vated by results from lVegetti & Koopmans! ( l2009l) . 
Equation l[T) can be easily extended to the case in which more than 
one lens is considered. Because the identification of a substructure 
in one system does not influence what we infer about another satel- 
lite in another lens, the Likelihood function for a set of ;i; lenses 
is simply the product of the independent likelihoods of individual 
detections. 



£{{n„m] I a,f,p) = ]~|i;(na.'nt I "J'P) 



(3) 



We now go through the details of the expectation value ii(a, f, < R), 
characterising the Poisson distribution of the number of substruc- 
tures in the potential of a generic lens galaxy. 



Statistics of mass substructure 3 



2.2 Substructure expectation value 

In the ideal case of infinite sensitivity to a mass perturbation, one 
would be able to recover the full mass function between Mmin and 
Mnax- In practice, only substructures with mass between M]„„ and 
Mhigh are observable, hence the expectation value for observable 
substructures is 



fi(a, f, <R)= fxoia, /, < R,p) 



dm 



(4) 



with the expectation from the full mass function given by 



tia{aJ,<R,p) : 



/(< R) Mdm(< R) 

Mmax 



r 



dm 



/(< R) Mdm(< R) 



(2-a) (mJ,7^ - AfjT^;') 
(l-a) (M^7^ - Afj;T^') 

tog(Mm„ / M„i„) 

tog(M„ax/Mmi„) 
(Mmax - Mmiii) 



or 2 , a 1 



a = 2 



(5) 



Hence, we assume the normalised true mass function to be given 
by a power-law 



dP 
dm 



(i-a) III ' 



log (Mmiix/A^min) 



ff # 1 



1 



(6) 



Mdm(< ^) and /(< i?) are the cumulative mass in dark matter and 
the cumulative fraction of dark matter in subhaloes within the con- 
sidered radius, respectively. However, the presence of noise on the 
data and the statistical uncertainty with which masses are measured 
introduce a scatter in the observed mass function, so that detections 
can be spread inside or outside our observational limits. The signif- 
icance of this effect depends on the substructure mass, with lower 
masses being affected by a larger relative uncertainty. The observed 
mass function dP/dm\co,n, can then be written as a convolution of the 
true mass function with the error distribution, which we assume to 
be Gaussian, hence 



lj(a,f,<R,p) = iio(a,f,< 



Jmi„„ dm 



dm — 



y"o(tt,/,< 



Jm,„„ JM„,i„ dm 



■ dm dm' . (7) 



A mathematical proof for this procedure can be found in the Ap- 
pendix. 



2.3 Posterior probability function of a and / 

Given a set of observations, in which a certain number of substruc- 
tures are identified and their masses are quantified for one or more 
lenses, equations l|T) and ^ can be used to infer the mass fraction 
and the mass function of the underlying subhalo population. Bayes' 
theorem relates the Likelihood function of the observations to the 
joint posterior probability of a and / in the following way 



P(aJ\{n„m],p): 



£({n,,m]\a,f,p)P(a,f\p) 
P({n,,m]\p) 



(8) 



where P{a,f\ p) is the prior probability density distribution func- 
tion of a and /. For the mass fraction we assume a non-informative 
uniform prior between the limits = and /n^ax = 1 ; while we 
test two different priors for a. We assume in one case a uniform dis- 
tribution between o-,,,,,, = 1.0 and a„„x = 3.0 and in the other case a 
Gaussian function with centre on 1.90 and a standard deviation of 
0.1, as found in almost all numerical simulations. We refer to the 
next section for a more detailed description of the prior probability 
density distributions. 

2.4 Dark matter mass 

As shown in equation l[5]l, the number of substructures expected 
in a given potential is a function of the dark mass Mdm(< R)- In 
this section we show how this mass can be empirically estimated. 
Specifically, we are interested in the cumulative mass within a nar- 
row annulus hR = 2 SR = 0.6" centred around the Einstein radius 
Re, where the formalism introduced can be considered valid. In the 
approximation of a small annulus, the dark matter mass contained 
in it can be approximated as 



Mdm ~ 4n Re SdmC^e) SR = An Re 



1 - 



2,0, SR . 



(9) 



E,o, and E, are the total and the stellar projected mass den- 
sity, respectively. It has been found by many authors, that the 
tota l mass density has a profile which is close to isothermal 



(e.g Gerhard et al.l200ll: Koopmans & Treull2002l : lKoopmans et"al 



l200dlCzoske et al.ll200alKoopmans et al.l2009l) . Similarlv. a llaffe 



19831) profile approximates the stellar mass distribution well, 

Ej. Re 



and 



2R 



(10) 



I.,(R) = C 



{i-^[(i-^-T'- (i-^T"(2-r>osh-(r') 



(11) 



with f = Rjr^. 

S(. is the critical surface density for tensing and r, is the scaling 
radius of the laffe profile, which relates to the effective radius of 
the lens galaxy as r, = RJOIA. We have assumed in the above 
equations that the Einstein radius and the effective radius are 
related to each other by R,, = 2Re, which is appr oximately the case 
for the average SLAGS lens jSolton et al.ll200M) . Obviously, this is 
not exactly true for any of these lenses, but we are not interested 
in analysing an exact reproduction of the SLAGS sample but just 
an average realisation of it. This can also be considered as a fair 
realisation of a typical massive early-type galaxy, given that both 
the internal and environmental properties of the SLAGS lenses do 
not significantly depart from those of other early-type galaxies 
with comparable velocity dispersio n and baryonic properties 
( lBoltonetai]|200d l Treu et al. I l2009l) . This assumption is in any 
case justified and does not influence our results; in fact, only the 
cumulative dark matter mass that is probed by all lenses is of 
relevance and its average value is not altered by the assumption of 
Rg = 2Re. 

The normalization constant C = 0.742^ can be derived by imposing 
E,o, — » S. a ^ for — > 0, that is by imposing that asymptotically 
for — > the mass density becomes that of stars only, i.e. the 
maximum bulge assumption. 

Because we assume that Re/r, is a constant and RJRe is also a con- 
stant, on average the projected DM mass fraction within an annulus 
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of 26R around the Einstein radius is a constant with a value of about 
63%. Hence we find that on average Edm « 0.5 x 0.63 at the 
Einstein radius for all lenses. T his fraction is consistent with other 
obser vations (e.g. iGavazzi eTal . 2008; Schechter & Wambsgans^ 
The number of substructures thus becomes only a function 
of the size of the aperture and the critical density, which itself only 
weakly depends on the source and lens redshifts. 



3 DATA REALISATION AND ANALYSIS 

In this section, we use a series of mock data sets to show that the 
formalism presented here allows us to constrain, in different situa- 
tions, the properties of the substructure population. Different data 
sets of lens galaxies with «/ = 10, 30 or 200 are analysed. A sam- 
ple based on the masses and radii of the SLACS lenses, ranked 
from the highest to the lowest mass enclosed within R^, is used 
to constru ct the dark matter mass func tion, while results from the 
method in IVegetti & Koopmans! l [2009h are used to set the observa- 
tional limits and the substructure mass uncertainty. In each lens, the 
substructures are distributed in mass according to dP/dm\„.ue with 
Mimn = 4.0 X 10'^ Mo and M,^^^ = 4.0 x 10' Mq ( Diem and et"al] 
12007 al ibi) and fluctuate in number with a Poisson probability distri- 
bution of expectation value y u(a, /, < R), given by equatio n lO. 
Although it was shown by IVegetti & Koopmana ( l2009h that the 
detection probability is a joint function of the substructure mass 
and position, this probability can be assumed to be independent 
of the perturber's position if a sufficiently small annulus around 
the Einstein ring is considered, in which case the equations pre- 
sented in the previous sections hold. We consider a region of ±0.3 
arcsec around Re, over which a typical SLACS lens shows a rea- 
sonable surface brightness of its images and within which we 
might expect to detect CD M substructures using the method of 
IVegetti & KoopmansI Jioogh . The extent of this area is, also, such 
that the mass fraction of substructures can be considered constant 
in radius, hence f(R) = const, over AR. 

Each mock data set is characterized by a different fraction /imc of 
substructures, while the true slope of the mass function is kept fixed 
at tttruc = 1 -90, as suggested by numerical simulations. Results from 
the latest high-resolution numerical simulations seem to indicate a 
dark matter mass fraction in substructures within a cylinder of pro- 
jected 10 kpc which is bet ween 0.3^ (Mark Vo gelsberger, private 
communication) and 0.5% jPiemand et al.l2008h . We therefore dis- 
cuss three different cases /true = 0. 1, 0.5, 2.5%. The latter high frac- 
tion is included because it is close to that suggested by the median 
value inferred from flux-ratio anomaly studies aPalal & Kochaneki 
I2OO2I) . 

3.1 Observational limits on the substructure mass 

We explore the effect of different values of the lowest detectable 
mass Miow, which is set at the statistical threshold above which we 
are confident that other effects do not create to o many false events. 
We have shown in lVegetti & KoopmansI ( l2009i) that, given the cur- 
rent HST data quality of the SLACS lenses, a lower mass limit for 

' This value has been obtained by including only particles within S2oo„„, 
adding masses of all those particles that are in subhaloes within the con- 
sidered cylinder and considering the median over 100 projected directions. 
No extrapolation beyond the resolution limit or cut on the subhalo mass is 
involved. 



a significant detection can be set around 10**Mo, depending on how 
close the perturbers are located with respect to the lensed images 
and the structure of the lensed images. However, these limits have 
been determined for cases not affected by systematic errors; we 
adopt Miow = 0.3, 1.0, 3.0 x 10** Mq. We set a finite upper mass fimit 
Mhigh = M„,,„. 

3.2 Priors on a and / 

While the mass fraction of satellites is the most uncertain param- 
eter, most studies seem to agree on the mass function, with val- 
ues of the slope a ranging from 1.8 to 2.0 (e.g.lHelmi et alj|2002l ; 
iGao et al.ll2004l ; lDiemand et al.ll2008l ; ISpringel et alj|2008l) . This is 
a direct consequence of the assumed cold nature of the dark mat- 
ter particles. We analyse two scenarios; the first, relying on results 
from numerical simulations, assumes for a a Gaussian prior centred 
at Q-truc with standard deviation it„ = 0.1; the second scenario, al- 
lowing for more freedom, has a uniform prior between 1.0 and 3.0. 
In general, the first case can be seen as a test of N-body simulations 
and the second as a test of nature itself, although in the specific 
case of this paper the data have been created with a combination 
of fraction and slope typical of a standard cosmology. The former 
prior, obviously, provides a well-defined mass-function slope, but 
it also reduces the uncertainty in the mass fraction, which is less 
well-defined than a and can vary considerably between similar sim- 
ulations. Different phenomena can affect the mass fraction, as for 
example the resolution of the simulations, or the lack, in high res- 
olution simulations, of gas physics, that could sensibly influence 
the substructure survival (i.e. / could be even higher than what the 
simulations suggest in the inner regions of the galactic haloes). We 
assume, conservatively, for / a uniform prior ranging between 0% 
and 100%. 



3.3 Results in the limit of no mass measurement errors 

We present results for cases in which the errors on the mass 
measurements can be neglected. Specifically, this translates into 
convolving the mass distribution not with a Gaussian but with a 
delta function around ni, in the equations of Section 2. 
In Fig.[T]we show the joint probability contours P(a,f\ {«,, m],p) 
and the marginalized probabilities P (f \ {n„m],p) and 
P (a \ {n,,m],p), for systems containing 10 randomly realised 
lenses. Specifically the plotted contours contain, in the limit of a 
Gaussian distribution, respectively, 68%, 95% and 97.2% of the 
marginalized probability function. 

In the case of a uniform prior, while a good upper limit 
to / can always be set, little can be said about the slope, 
which can only be constrained for a limited number of 
favourable physical and observational conditions, such as 
(/ = 0.5%, M,o„ = 0.3 X 10** Mo), (/ = 2.5%, = 0.3 x 10** Mq) 
and (/ = 2.5%,Mi„„ = 1.0 x 10** Mq). In Fig. |2] an equivalent 
plot is presented for systems with n; = 30; although even more 
stringent limits can be given for /, we are still unable to recover 
the underlying mass function for most of the possible scenarios. 
The situation can be substantially improved by increasing the num- 
ber of detectable substructures with a larger number of lenses. To 
provide insight into future capabilities, results from three samples 
of 200 lenses with / = 0.5% and respectively Miow = 0.3 x 10** Mq, 
Miow = 1.0 x IO^Mg, and Mio„ = 3.0 x 10'*Mo are given in 
Fig. [3] Currently, no uniform sample with 200 lenses with high 
signal-to-noise ratio and high resolution (equivalent to that of 
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HST) is available. However, forthcoming ground and space-based 
instruments (e.g LSST/JDEM, EVLA, e-Merlin, LOFAR and 
SKA) can provide these numbers (in fact beyond these) in the 
coming 5~10 years and the required data quality by dedicated 
follow-up. A detailed charac terisation of the CDM mass-f unction, 
through the technique of IVegetti & Koopmansl ( |2009|) , could 
therefore be realisable in the coming years if investments are made 
in large high-resolut ion and high-sensitivi ty lens surveys with 
these instruments (see lKoopmans et alj2009h . As can be seen from 
both Figs.[T]and[2] in the case of a Gaussian prior on a, tight limits 
can be set on the mass fraction for all possible combinations of the 
considered parameters. 

The results from this section are summarised in Table 1, where 
we report the input values for each parameter, the recovered 
maximum-posterior values (/mp;Q'mp) at which P(a,f \ jiij,»ij,p) 
reaches its maximum, the median, the 68% and 95% (ixsg and crgj), 
confidence levels of the marginalized probabilities P{f \ {/ij, m},p) 
and P (or I [ris, m},p) for both the cases of a uniform and a Gaussian 
prior on a. cr^g, and 0-95, in the particular case of a Gaussian 
distribution, respectively represent the Icr and 2cr error. 



3.4 The effect of mass measurement errors 

We explore now how the presence of uncertainty in the mass mea- 
surements affects our analysis. In particular, we consider three 
cases: cr,„ = 0.1 x 10* Mq with Miow = 0.3 x 10* Mq, cr,„ = 
1/3 X lO^Mo with Miow = lO^Mo, and cr„, = 1.0 x 10*Mo 
with Miow = 3.0 X 10* Mo i.e. a limiting signal-to-noise ratio of 
Mi(,„/cr„, = 3. The lens systems analysed here have n; = 10, 30, 200 
lenses and a mass fraction in substructures /ti„e = 0.5%. Relative 
likelihood contours are plotted in the three panels of Fig. O These 
have to be compared with the equivalent no-error results in Figs.[T] 
|2]and[3] Results are reported in Table 2. 

The effect of measurement errors on the substructure mass de- 
pends on the form of prior adopted for a, with uniform priors being 
more strongly affected than Gaussian ones. Errors as large as cr,„ = 
1/3 X 10* Mq combined with mass threshold of Mio„ = 10* Mq, can 
significantly influence even systems with 200 lenses . Specifically, 
these systems were created by drawing masses between M„ii„ and 
Mnax. then scattering each mass with a Gaussian distribution (i.e. 
mimicking the measurement uncertainty) and then using only those 
objects that fall within the detection range [Mio„, Mi,igi,j to constrain 
the fraction and the mass distribution. We refer to the Appendix for 
a mathematical proof that this way of proceeding is equivalent to 
drawing between Mio„ and Mhijh with a Poisson probability density 
distribution of expectation values given by equation l|4]l. 



4 CONCLUSIONS 

We have introduced a statistical formalism for the interpretation 
and the generalisation of subhalo detection in gravitational lens 
galaxies, that allows us to quantify the mass fraction and the mass 
function of CDM substructures. Given mock sets of lenses, with 
properties typical of a CDM cosmology, we have analysed how 
well the true parameters can be recovered. The formalism de- 
pends on several parameters, such as e.g. the number of lenses, 
the mass detection threshold and the measurement errors. It has a 
very general nature and, in principle, it could be used to statistically 
analyse substructure detection by flux ratio anomalies or timede- 
lay/astrometric perturbations as well. In practice, these methods 



would first need to show that their mass estimates are meaning- 
ful and second they would have to determine the probability distri- 
bution of flux-ratio anomalies or perturbations, either as function 
of the lens geometry or marginalized over all model parameters, 
which could be rather computationally expensive. The method has 
been tested on several mock data sets, with parameter settings based 
on our knowledge of the SLACS lenses. Several physical and ob- 
servational scenarios have been considered. We list here the main 
results: 

• If the number of arc/ring lens systems is <c 100, as is the 
case for current surveys (e.g. SLACS), the ability to constrain the 
mass fraction and the mass function of satellites still depends on 
the form of prior which is assumed for a. In particular, if results 
from numerical simulations are assumed to hold and a Gaussian 
prior with a^ac = 1 .9 ± 0. 1 is adopted, we are able to constrain both 
a and / for any data sets containing a number of lenses > 10, 
with improved limits for either increasing mass fractions, decreas- 
ing detection threshold or increasing number of lenses. If instead 
a wider range of possibilities is explored by assuming a uniform 
prior, one can still set strong limits on /, even for values as low 
as / = 0.1% and a detection threshold Mio„ = 0.3 x 10* Mq, but 
the mass function slope can be recovered only in a limited number 
of favourable cases, characterized by high mass fraction and low 
detection threshold. 

• Our ability to constrain a could be considerably improved ei- 
ther by increasing our sensitivity to substructures, i.e. by increasing 
the quality of the data, or by increasing the number of analysed ob- 
jects. Although competing with the quality of HST seems at the 
moment difficult, future surveys such as LSST/JDEM in the optical 
and EVLA, e-MERLIN, LOFAR and SKA in the radio, will surely 
lead to an increase in the number of known lenses by several or ders 
of magnitude (see lKoopmans et alj2009l : lMarshall et alj2009l) and 
dedicated optical and/or radio follow-up could provide equivalent 
or better data quality than HST. We expect therefore, in the foresee- 
able future to be able to characterise the galactic subhalo population 
with stringent constraint, both on the mass fraction and slope. 

• Although we have not explicitly performed a model compar- 
ison between different cosmologies, as for example CDM versus 
Warm Dark Matter (this would require an extra marginalization 
of the parameter space), the formalism introduced here, combined 
with the sensitivity of our method to CDM substructures, will al- 
lows us, in the future, to discriminate among these two scenarios 
and thus test the physics of dark matter. 
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Figure 2. Similar to Fig.[T]for systems witli 30 lenses. 
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Figure 3. Results for three samples with 200 randomly realised lenses with / = 0.5% and respectively M[„„ = 0.3 X 10** M© (left panel), Miow = 1 -0 X 10^ Mq 
(middle panel) and Mio„ = 3.0 X 10** Mq (right panel). The joint probability P(a, f \ {«,, in],p) contours and marginalized probabilities /"(/ | {/jj, m],p) and 
P(a I {«,, m],p) for a uniform prior (solid lines) and for a Gaussian prior in a (dashed lines) are shown. 
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Table 1. Results for systems with no mass measurement error for a uniform (left) and a Gaussian (right) prior on a. The input parameters are given in columns 
(1) to (3): the true mass fraction in substructure, /i,ue. the lower detection mass threshold, Mi„„, and the number of lens systems in a given sample, n;. The 
maximum posterior values of / and a are given in columns (4) and (5). Columns from (6) to (8) and (9) to (1 1) list the median, the 68% and 95% CL for / 
and a respectively. 



/true 

(%) 


(lO'^Mo) 


ni 


/mp 
(%) 


«MP 


(%) 


0"/,68 
(%) 


0"/,95 
(%) 


'^med 




68 




95 


0.1 


0.3 


10 


0.1410.14 


1.91 1 1.90 


0.18 1 0.15 


-0.06 1 -0.04 
+0.11 1 +0.06 


-0.09 1 -0.07 
+0.22 1 +0.10 


1.88 1 1.90 


-0.37 1 
+0.48 1 


-0.05 
+0.05 


-0.58 
+0.81 


-0.09 
+0.09 






30 


0.11 1 0.10 


1.691 1.90 


0.13 1 0.10 


-0.04 1 -0.02 
+0.08 1 +0.03 


-0.07 1 -0.04 
+0.15 1 +0.05 


1.641 1.90 


-0.25 1 
+0.29 1 


-0.03 
+0.03 


-0.40 
+0.51 


-0.06 
+0.05 



0.5 



1.0 



3.0 



0.3 



1.0 



3.0 



2.5 



0.3 



1.0 



3.0 



10 
30 

10 

30 



10 
30 

10 
30 



0.55 
0.10 

1.53 
0.80 



10 0.45 

30 0.49 

200 0.51 

10 0.52 

30 0.42 

200 0.46 

10 1.85 

30 0.42 

200 0.50 

10 2.50 

30 2.39 



2.42 
2.12 

2.74 



0.16 
0.12 

0.10 
0.13 



2.46 
2.14 



2.70 
1.81 

2.77 
2.62 
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0.50 1.92 

0.48 2.04 
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0.48 1.82 

0.46 2.56 

0.50 1.68 

0.50 1.90 

2.56 1.88 

2.31 1.82 



1.75 
1.70 



2.76 1. 



1.8612.48 1.2811 



1.90 
1.90 

1.90 
1.90 



0.60 
0.20 

2.15 
1.58 



1.89 0.49 

1.90 0.51 
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0.13 
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0.42 
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0.52 
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0.44 
0.50 

0.60 
0.44 
0.43 
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2.33 
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1.9612.53 
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-0.10 
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+ 1.40 


-0.10 
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+0.27 


-0.12 
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-0.04 
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+0.03 


-0.09 
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-0.01 
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+0.05 
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APPENDIX A: 

We show here that the procedure of first drawing objects between 
Mmin and M^ax from dP/dm\,n,e, then scattering with a Gaussian 
and finally restricting to those masses between Mio„ and Mhigh gives 



a probability P(A'obs) of observing A'obs objects that is equivalent to 
the Poisson probability density distribution, with expectation value 
/jj, of Nabs objects between Mio„ and M^,^ expressed by the convo- 
lution in equation l|7j. 

Let us divide the mass ranges [Mmin, Mm^x] and [Miow, Mhigh] re- 
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Table 2. Results for systems with a mass measurement error of 0.1 X 10^ A/©, 0.3 X 10* Mq, and l.Ox IO^A/q for a uniform (left) and a Gaussian (right) prior on 
a. The input parameters are given in columns (1) to (3): the true mass fraction in substructure, /truc> the lower detection mass threshold, Miow and the number 
of lens systems in a given sample, n;. The maximum posterior values of / and a are given in columns (4) and (5). Columns from (6) to (8) and (9) to (1 1) list 
the median, the 68% and 95% CL for / and a respectively. 



/true 'Wlow "/ /mP ttMP /mcd 0"/,68 a"/,95 "med 0"ff,68 0"ff.t)5 

(%) (l0'*Mo) (%) (%) (%) (%) 



0.3 


10 


0.53 


0.61 


2.24 


1.91 


0.57 


0.64 


-0.09 
+0.09 


1 -0.12 
1 +0.09 


-0.15 
+0. 18 


-0.15 
+0. 18 


2.25 


1.91 


-0.20 
+0.22 


-0.02 
+0.04 


-0.32 
+0.36 


-0.04 
+0.06 




30 


0.53 


0.53 


1 m 
1.91 


1.90 


0.54 


0.54 


-0.09 
+0.09 


1 -0.06 
1 +0.06 


-0.12 
+0.15 


-0.09 
+0.09 


1 fii 
1.91 


1.91 


-0.12 
+0.12 


-0.01 
+0.01 


-0.18 
+0.20 


-0.04 
+0.01 




200 


0.51 


0.52 


1.93 


1.90 


0.51 


0.51 


-0.06 
+0.06 


1 -0.03 
1 +0.06 


-0.09 
+0.12 


-0.06 
+0.09 


1.93 


1.91 


-0.08 
+0.08 


-0.01 
+0.01 


-0.14 
+0.14 


-0.02 
+0.01 


1.0 


10 


0.53 


0.44 


2.18 


1.90 


0.82 


0.48 


-0.36 
+0.94 


1 -0.15 
1 +0.15 


-0.48 
+ 1.79 


-0.21 
+0.27 


2.41 


1.91 


-0.46 
+0.38 


-0.04 
+0.04 


-0.77 
+0.52 


-0.08 
+0.06 




30 


0.55 


0.54 


1.98 


1.90 


0.64 


0.54 


-0.15 
+0.18 


1 -0.09 
1 +0.12 


-0.21 
+048 


-0.15 
+0.21 


2.03 


1.91 


-0.24 
+0.26 


-0.04 
+0.02 


-0.40 
+0.44 


-0.04 
+0.04 




200 


0.67 


0.59 


2.12 


1.90 


0.72 


0.61 


-0.15 
+0.18 


1 -0.09 
1 +0.09 


-0.21 
+0.39 


-0.12 
+0.15 


2.15 


1.91 


-0.16 
+0.18 


-0.01 
+0.01 


-0.28 
+0.30 


-0.04 
+0.02 


3.0 


10 


4.69 


0.42 


2.85 


1.90 


4.01 


0.50 


-2.76 
+4.76 


1 -0.25 
1 +0.25 


-3.51 
+9.02 


-0.25 
+0.50 


2.78 


1.91 


-0.38 
+0.16 


-0.05 
+0.06 


-0.83 
+0.20 


-0.09 
+0.09 




30 


2.76 


0.51 


2.57 


1.90 


5.26 


0.75 


-3.51 1 
+4.76 1 


-0.25 
+0.1e-4 


-4.26 
+8.02 


-0.25 
+0.25 


2.76 


1.91 


-0.32 
+0.17 


-0.04 
+0.03 


-0.62 
+0.21 


-0.06 
+0.06 




200 


0.51 


0.51 


1.90 


1.90 


3.76 


0.75 


-2.25 1 
+4.26 1 


-0.25 
+0.1 c-4 


-3.01 
+7.02 


-0.25 
+0.25 


2.68 


1.90 


-0.33 
+0.22 


-0.03 
+0.03 


-0.58 
+0.28 


-0.04 
+0.05 



spectively in n and n' sub-intervals of infinitesimally small widths 
dm and dm' . Lets call p^j the probability that an object is scattered 
from the j-\h mass bin dmj to the fc-th one dm^.. In the particular 
case of Gaussian errors, p^j^k reads as follows 

g-(m-m')"/2o-- 

Ps,i.k = = dm[ , (Al) 

V2ncr 

with m 6 dmj and m' e dm'i^. 

First we show that, if substructures are Poisson distributed in each 
mass bin with expectation value dp j, then the probability of having 
A', objects scattered from dm^ into rfni^ is also a Poisson distribution 
with expectation value dp^jj, = Ps.j,k dp j 

p,,j(N„dp,,j,k) = J p(i,dpj) U pj;;, (i - p,,,,,)""' = 



i=JVs 



if 1™Z 



dp) 



, 1 - Ps,j.k I 



J.N, ,,-N,, .l,,k 

dp. dp. k 

i=0 



{p.,.,kdpif 



(A2) 



hence it follows that dp^ jk = Psjj, dpj. This result directly follows 
from the fact that in the case of high number statistics the Binomial 
tends to a Poisson distribution and from the product rule of the 
Poisson distribution. Specifically, each dpj reads as 

dP 



dpj = po{a,f,R) 



dm 



dmi 



(A3) 



We now extend this result to two mass intervals of the same size 
dm; thanks to the sum rule, the probability of A'o^., objects being 



scattered is again Poissonian with dps^k = (Ps,i,kdpi + Psxkdpi) 

P(N„,„,dp,) = ^ P.,A')Ps.2iNob.. - = 
1=0 



Nob., 



1 1 ^ Mi*, V 

Nobs 



N„l„ 



-d^.,.t dp'^ob' 



(A4) 



By induction it can be shown that in the case of a generic number n 
of intervals dp^j^ = 2','=i P.!.j.k dpj, so that the probability of being 
scattered outside [M,,,;,,, M,^i,x] and inside [Mio„, Mhigh] is a Poisson 
distribution with expectation value 



k=l j=l k=l j=l 



dmj dm'f. . 



(A5) 



In the limit of equally infinitesimal intervals, i.e dm — > and 
dm' — » and making use of equation l lAIb this becomes 



q.e.d 



dm 



dm dm' . 



(A6) 



© 2002 RAS, MNRAS 000,[T]-?? 



